% Compute the global solution using a paramerized expectations algorithm
% Dynamic bank run model
% Bank leverage choice only
% Capital depreciation shock

close all
warning off

load('no_policy')

%% Generate random shocks

options = optimset;
T = 100000;               % simulation periods
t_drop = 1000;

rng('default')
rng(1)                   % same random shock 
shocks = normrnd(0,1,T+1,1);  % random shock
%Exi = zeros(T+1);  % stochastic steady state

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);
erk  = zeros(T+1,1);
x    = zeros(T+1,npol);
x_nonlog    = zeros(T+1,npol);
x_allnonlog = zeros(T+1,npol);

%% Simulation

zthre   = 0;

kv(1)  = stss_k;
rbv(1) = stss_rb;
nv(1)  = stss_n;
lv(1)  = kv(1)/nv(1);
av(1)  = 0;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));

        K_nonlog  = [log(kv(t-1)) rbv(t-1) log(nv(t-1)) xiv(t)];
        x0_nonlog = (K_nonlog).^expmx;
        x_nonlog(t,:)  = prod(x0_nonlog,2)';
        
        K_allnonlog  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0_allnonlog = (K_allnonlog).^expmx;
        x_allnonlog(t,:)  = prod(x0_allnonlog,2)';

        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x(t,:)  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);

            if pibv(t) > 0
                rhseev(t) = exp(x(t,:)*solved_eta1_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta1_rbar);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x(t,:)*solved_eta3_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta3_rbar);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end

        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x(t,:)*solved_eta2_rhsee);
            rbv(t)    = exp(x(t,:)*solved_eta2_rbar);
            runv(t) = 1;
        end

        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);

        erk(t) = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
               + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

        rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk(t),lambda,gamma,delta,sigma_rk);   % deposit interest rate

        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp     = (rkhats_next - erk(t)) / sigma_rk;
        cdf = normcdf(z_temp);
        prv(t) = cdf;
        
    end
    % End of T period simulation

%% create regression variables

years        = 3;
reg_num      = 0;
reg_num_posi = 0;

lhs      = zeros(10,1);
lhs_prob = zeros(10,1);
depo_g1  = zeros(10,1);
depo_g2  = zeros(10,1);
depo_g3  = zeros(10,1);
depo_g4  = zeros(10,1);
depo_g5  = zeros(10,1);
depo_g   = zeros(10,1);
dy_g     = zeros(10,1);
depo_g_years_mean = zeros(10,1);
depo_g_years_root = zeros(10,1);

lhs_posi      = zeros(10,1);
lhs_posi_prob = zeros(10,1);
depo_g_posi   = zeros(10,1);
dy_g_posi     = zeros(10,1);

    for t=2:T+1
        
        if t>499 && t<T+1-4 && floor(t/4)==t/4 && sum(runv(t-(4*years):t-1))==0
%        if t>499 %&& sum(runv(t-(4*years):t-1))==0
            
            reg_num = reg_num + 1;
            
            if xiv(t)>=1 && xiv(t+1)>=1 && xiv(t+2)>=1 && xiv(t+3)>=1 %|| sum(runv(t-(4*years):t-1))>0
                lhs(reg_num) = 0;
            else
                lhs(reg_num) = 1;
            end
            
            lhs_prob(reg_num) = prv(t);
            
            for y=1:years
%                depo_g(reg_num,y) = (kv(t-1-((y-1)*4))-nv(t-1-((y-1)*4)))/(kv(t-5-((y-1)*4))-nv(t-5-((y-1)*4))) - 1;
                depo_g(reg_num,y) = (kv(t-1-((y-1)*4)))/(kv(t-5-((y-1)*4))) - 1;
                dy_g(reg_num,y)   = (kv(t-1-((y-1)*4))-nv(t-1-((y-1)*4)))/yv(t-1-((y-1)*4)) - (kv(t-5-((y-1)*4))-nv(t-5-((y-1)*4)))/yv(t-5-((y-1)*4));
            end
%            depo_g_years_mean(reg_num,1) = ( (kv(t-1)-nv(t-1)) / (kv(t-1-(years*4))-nv(t-1-(years*4))) - 1 ) / years;
%            depo_g_years_root(reg_num,1) = nthroot( (kv(t-1)-nv(t-1)) / (kv(t-1-(years*4))) ,years) - 1;
            depo_g_years_mean(reg_num,1) = ( (kv(t-1)) / (kv(t-1-(years*4))) - 1 ) / years;
            depo_g_years_root(reg_num,1) = nthroot( (kv(t-1)) / (kv(t-1-(years*4))) ,years) - 1;
            
        end
        
        % check if negative profit matters for the result
        
        if pibv(t-1)>0 && pibv(t-2)>0 && pibv(t-3)>0 && pibv(t-4)>0
        if t>499 && t<T+1-4 && floor(t/4)==t/4 && sum(runv(t-(4*years):t-1))==0
%        if t>499 %&& sum(runv(t-(4*years):t-1))==0
            
            reg_num_posi = reg_num_posi + 1;
            
            if xiv(t)>=1 && xiv(t+1)>=1 && xiv(t+2)>=1 && xiv(t+3)>=1 %|| sum(runv(t-(4*years):t-1))>0
                lhs_posi(reg_num_posi) = 0;
            else
                lhs_posi(reg_num_posi) = 1;
            end
            
            lhs_posi_prob(reg_num_posi) = prv(t);
            
%            depo_g1(reg_num) = (kv(t-1)-nv(t-1))/(kv(t- 5)-nv(t- 5)) - 1;
%            depo_g2(reg_num) = (kv(t-1)-nv(t-1))/(kv(t- 9)-nv(t- 9)) - 1;
%            depo_g3(reg_num) = (kv(t-1)-nv(t-1))/(kv(t-13)-nv(t-13)) - 1;
%            depo_g4(reg_num) = (kv(t-1)-nv(t-1))/(kv(t-17)-nv(t-17)) - 1;
%            depo_g5(reg_num) = (kv(t-1)-nv(t-1))/(kv(t-21)-nv(t-21)) - 1;
            for y=1:years
                depo_g_posi(reg_num_posi,y) = (kv(t-1-((y-1)*4))-nv(t-1-((y-1)*4)))/(kv(t-5-((y-1)*4))-nv(t-5-((y-1)*4))) - 1;
                dy_g_posi(reg_num_posi,y)   = (kv(t-1-((y-1)*4))-nv(t-1-((y-1)*4)))/yv(t-1-((y-1)*4)) - (kv(t-5-((y-1)*4))-nv(t-5-((y-1)*4)))/yv(t-5-((y-1)*4));
%                depo_g(reg_num) = (kv(t- 5)-nv(t- 5))/(kv(t- 9)-nv(t- 9)) - 1;
%                depo_g(reg_num) = (kv(t- 9)-nv(t- 9))/(kv(t-13)-nv(t-13)) - 1;
%                depo_g(reg_num) = (kv(t-13)-nv(t-13))/(kv(t-17)-nv(t-17)) - 1;
%                depo_g(reg_num) = (kv(t-17)-nv(t-17))/(kv(t-21)-nv(t-21)) - 1;
            end
%            depo_g1(reg_num) = (kv(t)-nv(t))/(kv(t- 4)-nv(t- 4)) - 1;
%            depo_g2(reg_num) = (kv(t)-nv(t))/(kv(t- 8)-nv(t- 8)) - 1;
%            depo_g3(reg_num) = (kv(t)-nv(t))/(kv(t-12)-nv(t-12)) - 1;
%            depo_g4(reg_num) = (kv(t)-nv(t))/(kv(t-16)-nv(t-16)) - 1;
%            depo_g5(reg_num) = (kv(t)-nv(t))/(kv(t-20)-nv(t-20)) - 1;
            
        end
        end

    end

%% OLS

%rhs         = horzcat(depo_g1,depo_g2,depo_g3,depo_g4,depo_g5);
%rhs_sum     = depo_g1+depo_g2+depo_g3+depo_g4+depo_g5;

rhs     = depo_g;
rhs_sum = sum(depo_g,2);

reg_results      = regstats(lhs,rhs);
reg_results_sum  = regstats(lhs,rhs_sum);
reg_results_prob = regstats(lhs_prob,rhs);

reg_results.beta

stdv_depog = std(depo_g(:,1));
stdv_depog_years_mean = std(depo_g_years_mean);
stdv_depog_years_root = std(depo_g_years_root);
runprob_effect = stdv_depog_years_root * ( reg_results.beta(2)+reg_results.beta(3)+reg_results.beta(4) )

%% Logit

[reg_logit_beta,reg_logit_dev,reg_logit_stats] = glmfit(rhs,lhs,'binomial','link','logit');

reg_logit_beta

mean_prob = mean(lhs);
logit_eval = zeros(1,3);
logit_eval(1) = reg_logit_beta(2) * mean_prob * (1-mean_prob);
logit_eval(2) = reg_logit_beta(3) * mean_prob * (1-mean_prob);
logit_eval(3) = reg_logit_beta(4) * mean_prob * (1-mean_prob);
logit_eval
